The most important DeepBlue entity is the experiment. The experiments contains the Epigenetic Data with their data description, that we name as metadata.
The most important DeepBlue entity is the experiment. The experiments contains the Epigenetic Data with their data description, that we name as metadata.
We can list the experiments with the list_experiments command. The command has 6 parameters: genome, epigenetic_mark, sample, technique, project, user_key, which only user_key is mandatory and the others parameters work as filters.
import xmlrpclib
user_key = raw_input("Inform your user_key")
url = "http://deepblue.mpi-inf.mpg.de/xmlrpc"
server = xmlrpclib.Server(url, encoding='UTF-8', allow_none=True)
print server.echo(user_key)
From ours previous examples, lets list all experiments with the Epigenetic Mark H3K27ac of the ENCODE and BLUEPRINT projects:
(status, experiments) = server.list_experiments("", "H3K27ac", "", "", ["ENCODE", "Blueprint Epigenetics"], user_key)
We can list the experiments that contains Samples from blood:
(status, samples) = server.list_samples("blood", {}, user_key)
len(samples)
Obtain the sample IDs:
samples_ids = [sample[0] for sample in samples]
(status, experiments) = server.list_experiments("", "H3K27ac", samples_ids, "", ["ENCODE", "Blueprint Epigenetics"], user_key)
experiments
In this case, we are only getting the experiments where the Sample has exactly the BioSource blood. We could use the get_biosource_children:
(status, blood_children) = server.get_biosource_children('blood', user_key)
biosource_names = [biosource[1] for biosource in blood_children]
(status, samples) = server.list_samples(biosource_names, {}, user_key)
len(samples)
samples_ids = [sample[0] for sample in samples]
(status, experiments) = server.list_experiments("", "H3K27ac", samples_ids, "", ["ENCODE", "Blueprint Epigenetics"], user_key)
experiments
len(experiments)
Now, lets asume that we want all experiments that contains peaks. For that we can use the info command to gather more information about the experiment. The info command can be used to obtain more information about the experiments:
server.info(peak_experiments[0][0], user_key) # ID of the first experiment
We will filter all experiments that contains peaks verifying if the metadata field "data_type" is "peak" and (The other possible value for this field is "signal".) and the experiment constains the columns: CHROMOSOME, START, END, NAME, SCORE, STRAND, SIGNAL_VALUE, P_VALUE, Q_VALUE, that is the BED 6+3 format.
peak_experiments = []
peak_format = "CHROMOSOME,START,END,NAME,SCORE,STRAND,SIGNAL_VALUE,P_VALUE,Q_VALUE"
for experiment in experiments:
(status, info) = server.info(experiment[0], user_key) # experiment[0] is the ID of the experiment.
if info[0]["data_type"] == "peaks" and info[0]["format"].startswith(peak_format): # We use startwith() because we also want experiments that contain more columns.
peak_experiments.append(experiment)
len(peak_experiments)
The first step is to make a list of all experiments names. (Using the IDs is also possible, but it is recomended to use the names).
experiments_names = [experiment[1] for experiment in peak_experiments]
We can select the regions, for that we must use the command select_regions that has the following parameters: experiment_name, genome, epigenetic_mark, sample_id, technique, project, chromosome, start, end, and user_key. As we already have the experiments names, we can just use this list. We also will retrieve only the regions from the chromosome 1 from the starting position 0 to 100000000.
(status, query_id) = server.select_regions(experiments_names, None, None, None, None, None, "chr1", 0, 100000000, user_key)
query_id
We can obtain the regions using the get_regions command. This command requires the query_id, that was returned by the select_regions and the content format, where the @NAME is the experiment name, @BIOSOURCE the BioSource and @SAMPLE_ID the Sample ID of the Region.
format = "CHROMOSOME, START, END, NAME, SCORE, @NAME, @BIOSOURCE, @SAMPLE_ID"
(status, request_id) = server.get_regions(query_id, format, user_key)
request_id
The command get_regions returns a request_id. It is due the fact that all the data is processed dinamically, that provides all format flexibility. Some requests processed almost immediately, where others may take minutes or even hours. It depends mainly on the amount of data that will be processed as well the server usage. The request status can be verified with the get_request_status:
(status, request_status) = server.get_request_status(request_id, user_key)
print request_status
The request data can be retrieved when the request status be done.
(status, data) = server.get_request_data(request_id, user_key)
print format
print data